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1.  INTRODUCTION 


The  Impact  of  a  long,  slender,  eroding  rod  at  high  speed  on  a  thick  semi-infinite  target  was 
Initially  formulated  by  Alekseevski  (1966)  and  Tate  (1967,  1969).  The  governing  equations,  using  the 
notation  of  Wright  and  Frank  (1988),  are: 

i^V-U.  (1) 

LU^-Yip,,  (2) 

1/2  p,  (U-V/  K  =  //2  p,  +  R.  (3) 

and 

>  =  V.  (4) 

where  U  is  the  speed  of  the  tear  of  the  penetrator.  L  is  the  instantaneous  penetrator  length,  V  is  the 
penetration  velocity.  P  is  the  depth  of  penetration,  p,  is  the  penetrator  density.  Y  is  the  penetrator  yield 
stress.  Pi  is  the  target  density,  and  R  is  the  target  resistance.  In  the  equations,  a  dotted  quantity 
represents  the  dme  derivative,  d/dt. 

Wright  and  Frank  (1988)  and  Frank  and  Zook  (1987)  discuss  these  equations  in  detail,  including 
the  assumptions  made  in  the  derivation  and  approximate  solutions.  Our  intent  is  to  analyze  the 
mathematical,  not  the  physical,  aspects  of  the  equations  (1)  through  (4).  Basically,  L,  U,  V,  and  P  are 
the  unknown,  dependent  variables,  t  (the  dme)  is  the  independent  variable,  and  p,,  Y,  and  P  are 
known  constants. 

2.  THE  SOLUTION 

An  exact  analytical  solution  of  equations  (1)  -  (4)  for  L,  U,  V.  and  P  is  now  obtained.  From 
equation  (3) 

V  -  (5) 

1-y 


where,  if  is  the  initial  (known)  penetrator  velocity.  v=VAJo.  u=U/Uo,  T=p/p,,  and 
I=2(R  -  Y)/(p,  Uo^).  The  minus  sign  is  chosen  for  the  radical  in  the  solution  of  the  quadratic 


1 


equation  (3),  to  guarantee  that  V  remains  less  than  U.  with  both  U  and  V  real.  Note  that  for  the  case 
R  >  Y  (£  >  0),  the  minimum  admissabie  value  of  v  is  zero,  corresponding  to  the  moment  tiiat 
penetration  ceases,  which  occurs  u  \x=  For  the  case  of  R  <  Y  (£  <  0),  the  minimal  admissabie 
value  of  u  is  'l-t/f,  in  order  to  keep  the  root  real  in  (S).  In  this  case  v  a  u.  corresponding  to  the 
situation  where  rod  erosion  ceases  and  rigid  body  penetration  commences. 

Next,  from  equation  (2).  with  K  =  Y/(q,  u/), 

L  it  ^  “K  (6) 


Differentiation  gives 


LU-t  uLs  0, 


and  the  solution  for  dL/dt.  eliminating  L.  is 


Alternately,  from  equations  (1)  ai>d  (5). 


(ay -/yiJTzd-Y)  ). 


Combining  these  two  expressions  to  eliminate  dL/dt  gives 


U  _  yub  b^yu^  ♦  E(1  -y) 
b  "  K(l-y)  ■  A(l-y) 


Straightforward  integration  yields 


//i(ri{uv/Y  *^yu^*Ti\~y)  } 


yu^ 


2Ki\-y)  2K{l-y) 


/y«^  +  E(1  -y)  ^  G  . 


(7) 


(8) 


1 

I 


■ 

I 
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where  O  is  i  constant  of  integratitm  which  results  horn  evaluation  of  the  Integral  at  ttie  onset  of 
penetration,  when  u«  1  andOaUg.  Note  that  d  which  equals  (lAlo)<fU/dtl  and  has  dimensions 
of  (lAl.  can  be  evaluated  tirom  equation  (6)  as  U,  =  -K  The  coastant  G  may  be  expressed  as 


C^N*lnM  +  lnU^. 


where 


Af-tv/Y  ]  ; 

w.  ♦E(i-y)  _  Y 
2Ar(l-Y)  2Ar(l-Y) 


By  substituting  the  constant  A  for  the  exponent,  £/(2KVy).  equation  (8)  may  be  expressed  in  the 
following  fonn: 


*^yu^  -t-  Efl  -yT 

6,  M 


Y“* 

“  2K(\-y) 


2Kil-y) 


/y  w*  + 


1-Y)  *N 


(9) 


By  introducing  the  following  transfonnation  variable  z, 


•u^/y  -y ■ 


(10) 


equation  (9),  under  the  transformation  (10),  yields 


(jM-D/J  4  53(1  _y)7(^-5>'^)cxp(flz  -  Clz)dz  *Fdt, 


(11) 
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where 


and 


a- - 1 - 

%K^/yi^/y  ^l) 

P(1-Y)(Vy 

— 


F»4U^My/ytxp  N-^  . 


The  use  of  die  transfonnadon  (10)  produces  a  diHerendal  equation  (1 1).  in  which  the  variables,  z 
and  t,  are  now  separate.  If  integrated  from  dme  0  to  some  finite  dme  t  in  the  penetiadon  process,  the 
limits  on  z  will  vary  from  its  inidal  value,  when  u  =  1.  of 


r.  •  (/T  ♦Vy  ♦X(1-Y)  )  . 

to  some  intermediate  value  z.  For  the  case  of  R  >  Y.  the  terminal  value  of  dme  at  which  the 
governing  equations  are  applicable  occurs  when  penetradon  ceases  at  v  =  0,  in  which  case  u  ^  V  £ 
and  the  terminal  value  of  z,  expressed  as  z^,  is  given  by 


♦  l)  • 

For  the  case  of  R  <  Y,  the  long  rod  penetradon  equations  are  only  valid  (without  modification)  to  the 
dme  at  which  the  penetrator  begins  rigid  body  penetration,  in  which  case  u  =  v  =  V-2^,  and  the 
terminal  value  of  z  becomes 


Note  that  z  is  always  positive,  since  £  is  positive  for  R  >  Y,  and  (-£)  is  also  positive  for  R  <  Y. 
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Equation  (1 1)  may  be  Amber  simplifled  by  letting 


4.t 


in  the  first  integral  over  z  and 


0.2 


in  the  second  integral  over  z.  Under  these  transfbnnations,  and  letting  E]s2/(A-fl)  and  E}=2/(A-1).  the 
integration  of  equation  (11)  reduces  to 

* 


•  I 

E(  1  -  Y ) £,  Jexp  -  ce*^)  dG  -  £  Jdr . 

•.  0 


(12) 


The  solution  is  now  reduced  to  a  straightforward  integration,  though  it  requires  evaluation  of  an 
exponential  integral  which,  in  theory,  is  a  ubulated  function  of  five  input  parameters.  Defining  the 
function.  W.  as 

W'(£,C.£.y,.yj)»£  Jexp(By‘-C>'"*)dy.  (13) 


the  solution  for  t  becomes  simply 


r -  [W(£,C.£,,(|» )  ♦  £(  1  - Y ) W(£,C,£j.e,,0)]/F . 


In  practice,  our  function  is  evaluated  by  expanding  the  exponential  function  in  equation  (13)  in  a 
power  series  and  integrating  term  by  term,  to  the  desired  degree  of  precision.  The  result  is  z  as  an 
implicit  function  of  the  time  variable,  t. 


The  number  of  power  series  tenns  required  for  convergence  of  our  W  function  varies  a  great  deal 
with  the  input  conditions  to  the  problem.  In  particular,  the  evaluation  of  our  W  function  by  way  of 
power  series  can  be  exacerbated  for  problems  where  the  penetrator  velocity  overwhelms  the  strengths 
of  the  rod  and  target  materials.  Fortunately,  problems  in  this  velocity  range  are  generally  beyond  the 
range  of  interest,  for  typical  long  rod  penetrator  impacts. 

Because  B  is  always  positive,  the  first  part  of  the  exponential  term  grows  with  z.  As  B  is  made 
parabolically  larger  by  increasing  the  penetrator  striking  velocity  U^,  more  terms  are  required  to  make 
the  power  series  converge.  Because  it  is  a  binomial  that  needs  to  be  exponentiated,  use  of  n  terms  in 
the  exponential  expansion  requires  that  n(n-l)/2  monomials  be  evaluated.  Similarly,  the  coefficient  for 
each  of  the  n  highest  order  motK>mials  requires  (2n)  operations  to  evaluate.  Thus,  the  computational 
effort  required  to  evaluate  our  W  function  varies  greatly  with  initial  conditions  to  the  problem. 

Typical  penetration  problems  involving  significant,  but  not  total,  penetrator  erosion  require  that 
10  to  20  exponential  terms  be  evaluated  in  order  to  keep  the  relative  error  of  the  time  variable  in  the 
fifth  decimal  place.  Such  calculations  require  mere  seconds  of  computation  on  a  PC.  As 
hypervelocity  conditions  are  approached,  the  number  of  exponential  terms  required  for  the  same 
convergence  epsilon  may  exceed  200  (recall  that  2(X)  exponential  terms  implies  200  x  199/2  =  19,900 
monomials),  requiring  several  minutes  on  a  PC.  Fortunately,  this  solution  technique  need  not  be 
pursued  for  problems  in  hypervelocity  since,  under  these  conditions,  the  long  rod  penetration  equations 
approach  the  standard  Bernoulli  flow  conditions,  which  may  be  readily  solved  by  hand. 

Having  t  as  a  function  of  z,  the  normalized  rod  speed,  u,  follows  from 
equation  (10)  as 


2\/y 


(14) 


Equation  (5)  may  then  be  employed  to  obtain  the  normalized  penetration  velocity,  v.  The  rate  of  rod 
erosion  comes  from  equation  (1),  in  the  fonn 

L  ^  U^(y  -u) . 


f) 


The  ^<eaetTator  length,  L,  may  be  obtained  in  the  following  fashion.  From  equations  (1)  and  (2),  one 
obtains 


L  _  (ii-v)6 
L  JC 


Substituting  for  v  and  u  give  foe  following: 


L_  yuk  _  fiVvw*  +  -Y) 

I*  ^(1-Y)'  ^(1-y) 


Note  the  identical  fonn  of  this  relatior  and  equation  (7).  As  a  result,  the  solution  lor,,..,  nearly 
identical  to  equation  (9),  including  the  definition  of  constants  A,  M,  and  N: 


{a^ ♦  ^YaUE(l-Y)  } 


M 


^  yu^ 

2K{l-y) 


2/f(l-Y) 


/y“* ♦  E(1  -y) 


Finally,  the  penetration,  P,  is  obtained  as  follows: 


f  I 

P»  jvdfU^jv{di/dz)dz. 

o  o 


(15) 


The  quantity  dt/dz  has  been  previously  obtained  in  equation  (1 1)  as 


dt  _  +  5^(1  -y)z^^-‘>y^)cxp(Bz  -  C/z) 

Iz’  *  F 


We  may  express  v  in  terms  of  z,  using  equations  (5)  and  (14).  as  follows: 


V 


1 

2^ 


(i-,/7)'/r 

(i-y) 


r(wV7) 
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These  substitutions  into  (IS)  produce  the  following  expression  for  P: 


m 


l-'/f  _  jA/I  j(A-J)/l  _  g*(U>/Y)(l-Y) 

L2v^Y(1-Y)  2v/y 


exp(Bz  -  Clz)dz. 


This  equation  is  similar  to  Equation  11.  in  that  it  may  be  solved  directly  for  penetration  P,  in  this  case, 
in  terms  of  serveral  W  functions. 

The  results  of  the  present  solution  have  been  compared  with  the  original  results  presented  by  Tate 
(1967).  in  which  he  numerically  integrated  the  pefKtration>time  history  of  a  duralumin  rod  strinking  a 
polythene  target,  for  two  different  target  resistance  values.  The  curves  resulting  from  the  present 
analytical  solution  achieve  a  direct  overiay  to  Tate’s  numerically  integrated  solution. 

Finally,  the  source  code  listing  for  the  penetration  equations  considered  in  this  report  is  given  in 
the  Appendix. 

3.  SPECIAL  CASE  SOLUTIONS 
Two  special  cases  are  considered: 

A)  p,  =  p,  =  p.  and  Y  =  R  »  o 

B)  p,  =  p,  =  p. 

SPECIAL  CASE  A;  p,  =  p,  =  p,  and  Y  =  R  «  o. 

Equation  (3)  reduces  to 


(u-v)*  -v^, 

with  the  non-trivial  solution  u  -  2v.  Differentiating  equation  (2),  and  combining  the  result  with 
equation  (1),  as  before,  yields: 
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where 


a  2ua 

1"“^’ 


Integrating  this  equation  and  evaluating  the  constant  of  integration  at  time  equd  0,  where  d  - 
which  has  the  value  d, «  (-If  U^/(4L«).  results  in 

Integrating  again  for  u  gives  the  result  in  terms  of  our  W  function  as 


expCJ£0  0.2. 

lU, 


l.tf). 


which  gives  u  implicitly  as  a  function  of  time,  t  As  mentioned  above,  u  ■  2v  apices,  and  thus 
determines  penetration  rate  v.  Also, 


L  • 


-V  u 

9 

-y~ 


and 


L  ■ 


-HU. 


follow  directly.  To  determine  penetration  P,  employ  the  tactic  of  equation  (IS),  namely 

I  9 

p.jvdtm  U,jv(dt/du)  du . 

9  I 

Penetration  rate,  v,  is  known  directly  in  terms  of  u,  equal  to  (u/2),  and  dt/du  is  simply  1/u,  given  by 


dt  1 

Tu’rr\-f, 


1 


Tu^\nr 

Thus,  making  use  of  the  term  u.  ■  (-H  UJ/(4L..),  the  penetration  may  be  computed  as 


P^L. 
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SI^CIAL  CASE  B:  p,  =  p,  s  p. 


Equation  (3),  which  is  no  longer  quadratic,  as  In  tlu:  general  case,  yields 

V  ■  l(tt-IVM). 

2 

Differentiating  equation  (2),  and  combining  the  result  with  equation  (1),  and  the  expression  for  v 
above,  yields 

a  ^  _  uU  _  Ei 

l"  “  Tu' 

Direa  integration  results  in 

where  the  exponent  A  has  the  value,  analogous  to  the  general  case  derivation,  of  £/(2K).  The 
variables  are  separable,  and 


This  integral  is  evaluated  using  the  same  procedure  used  for  equation  (11).  By  letting 
the  integral  reduce  to 


(A 


L,fexp[il-.^J 


the  evaluation  of  which  may  be  expressed  in  terms  of  our  W  function  as 

r ' 


*  exp  w((4Ar)  ‘.0.2/(/lH).l.u<^*'>). 
26,  [4  A 


The  evaluation  of  the  remaining  variables  v,  dL/dt,  L,  P  follows  an  approach  analogous  to  the 
general  case. 
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4.  EXTENDING  THE  SOLUTION 


The  long  rod  equations  (l)-(4)  do  not  hold  for  the  complete  penetration  process.  In  paiticular,  if 
R  >  Y,  the  equations  are  only  valid  until  penetration  ceases  (V  =  0).  even  though  rod  erosion  is  still 
occuiing  (U  >  0).  On  the  other  hand,  for  R  <  Y.  the  equations  (l)-(4)  are  only  valid  until  rod  erosion 
ceases  (V  =  U),  even  though  rigid  body  penetration  proceeds  (V  =  U  >  0). 

When  either  of  these  limiting  conditions  occurs,  the  equations  (1)*(4)  can  no  longer  remain  valid 
without  some  sort  of  alteration,  in  order  to  prevent  the  situations  V  <  0  or  U  <  V.  To  illustrate  this 
point,  consider  adding  the  constraint  of  setting  V  to  zero,  when  R  >  Y,  to  obtain 

L  =-I/. 

L  0  =-r/p,. 

112  p,U^  +  Y  =  R,  and 
P  =  0 

or  setting  the  term  V  =  U,  if  R  <  Y,  to  obtain 

L  ^  0, 

LU  =  -r/p,, 

Y  *  112  p,  H  .»?.  and 

P*  U. 

Since  V  has  been  eliminated  as  a  dependent  variable,  leaving  just  U,  L.  and  P  as  unknowns,  the 
stipulation  of  four  governing  equations  overconstrains  the  problem.  If  no  other  conditions  are 
changed,  then  the  third  equation  of  each  set  (modified  Bernoulli  equation)  amounts  to  setting  the  rod 
and/or  penetration  velocity  to  a  constant,  which  is  clearly  incorrect.  One  possibility  involves  the 
elimination  of  the  modified  Bernoulli  equation  directly,  on  the  assumption  that  the  U,  V  relationship  it 
defines  is  replaced  by  the  rigidity  constraints:  V  =  0  if  R  >  Y,  or  V  =  U  if  R  <  Y. 

Another  possible  remedy  to  this  problem,  and  the  alternative  being  proposed,  involves  introducing 
an  additional  unknown,  in  the  form  of  a  material  resistance  R  or  Y,  to  vary  from  its  initial  value  of  R„ 
or  Yp,  respectively. 

The  rationale  for  this  step  is  given  now.  The  term  R  and  Y  represent  dcviatoric  stres.s  rc.sisianccs 
built  up  in  the  target  and  penetraior,  respectively,  which  act  as  part  of  a  force  balance  in  the  modified 
Bernoulli  in  equation  (3).  For  the  case  of  R^  >  Y^,  the  target  penetration  ceases  at  some  point  while 
rod  erosion  continues.  After  this  point  of  target  rigidity,  the  target  is  elastic.  Thus,  the  resistive  stress 


II 


offered  by  the  target  (R)  will  decrease  from  its  hill  i^astlc  value  (Rg).  to  an  elastic  value  exactly  equal 
to  the  penetiator  yield  strength  (Y^.  just  as  the  penetrator  velocity  becomes  identically  zero.  For  the 
alternative  case  of  R,  <  Y,,  the  rod  erosion  ceases  at  some  point  while  rigid  body  penetradon 
continues.  After  this  point  of  penetrator  rigidity,  the  penetiator  is  elastic.  Thus,  the  resisdve  stress 
offered  by  the  poietrator  (Y)  will  decrease  firom  its  full  plastic  value  (Y^)  to  an  elastic  value  exactly 
equal  to  the  target  yield  strength  OU,  just  as  the  penetration  velocity  becomes  identically  »ro. 

These  equation  seu  may  then  be  solved  to  detenrine  the  residual  rod  erosion  (if  R  >  Y)  or 
residual  pmietration  (if  R  <  Y),  and  the  associated  event  duration,  using  the  ai^prlate  initial 
conditicms,  wnich  resulted  from  the  terminal  conditions  associated  with  the  previous  solution  of 
Equations  (1)*(4).  The  boundary  coivlitions  are  given  below,  with  subscript  r  referring  to  conditions  at 
the  onset  of  rigidity  (Uiget  rigidity  if  R  >  Y  or  rod  rigidity  if  R  <  Y),  and  subscript  x  referring  to 
conditions  of  U  a  V  s  0. 


Target  Rigid:  R  >  Y  (V  ■  0) 


Initial  Conditions 

U  a  U, 

L  =  4 

P  a  P, 

t  =  t, 

R  =  Ro 

Y  =  Y, 


Final  C!ondilions 

U  =  0 
L  =  L, 

P  =  P, 
t  =  ^ 

R  =  Y, 

Y  =  Y, 


Penetrator  Rigid:  R  <  Y  (V  = 


Initial  Conditions 


a./  ••  Wj 

L  »  L, 

P  a  P, 

t  a  1, 

R  *  Ro 

Y  =  Y, 


Final  Conditions 


U  =  0 
L  =  L, 
P  a  P, 

‘  =  ^ 
R  =  Ro 
Y  =  Ro 


For  the  case  of  R  >  Y,  the  solution  is  expressable  in  terms  of  the  W  function  described  in  the 
original  section  of  this  report.  Non-dimensional  ternis  (e  g..  K,  y,  u)  defined  earlier  arc  also 
employed,  to  give  the  solution  as: 
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The  tenm  refere  to  the  value  of  K  at  its  original  constant  value,  when  R  =  R,  and  Y  =  Y„.  For  the 
case  of  R  <  Y,  the  solution  is  deteimined  as: 


2Y. 


|w 


yata 


tan 


-1 


U. 


T?r 


The  results  are  again  compared  to  Tate  (1967),  in  which  the  penetration  of  a  duralumin  rod  into 
polythene  was  modeled  for  two  values  of  target  resistance — 14  and  27  tons  per  square  inch.  Tate 
pointed  out  that  these  two  values  produced  curves  vdiich  straddled  the  actual  penetration*time  data.  In 
Tate’s  fomtulation,  the  equations  are  only  valid  to  the  point  where  rigid-body  penetration  commences 
which,  for  the  case  in  question,  accounts  for  only  80%  of  the  total  penetration.  Using  the  current 
formulation,  for  computing  the  residual,  rigid-body  penetration  of  the  duralumin  rod,  the  experimental 
value  of  penetration  may  be  used  to  compute  die  exact  value  of  target  resistance  required.  For  this 
case,  using  the  units  of  Tate,  a  target  resistance  of  20  tons  per  square  inch  was  computed  as  necessary 
to  produce  the  experimentally  observed  residual  penetration,  measured  from  Tate’s  graph  as  6.15 
inches.  The  penetration  process  is  predicted  to  stop  at  242  psec  alter  impact,  thus  implying  that  rigid 
body  penetration  takes  42%  of  the  toul  time  of  penetration. 


5.  CONCLUSIONS 


An  analytical  solution  to  the  long  rod  peneuation  equations  for  long  tod  penetration  is  offered. 
The  general  case  is  solved  as  well  as  two  special  cases  in  which  some  of  the  target  and  penetrator 
parameters  (e.g.,  density  and/or  strength)  are  equal.  This  analytical  solution  allows  a  faster  and  easier 
solution  of  the  penetration  equations,  since  stability  considerations  associated  with  any  numerically 
integrated  solution  ate  avoided. 

Additionally,  the  proposed  modification  to  the  original  model  permits  the  computation  of  residual 
rod  erosion  and  residual  penetration,  which  offer  additional  information  about  the  penetration  process, 
not  available  in  the  original  penetration  equations  formulated  by  Tate. 
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program  tata 

c 

C  CODED  BY  staven  B.  Sagletas,  AUGUST-SBPTBMBER,  1990. 

C 

C  THIS  PROGRAM  IMPLEMENTS  WALTERS'  MID  SE6LETES'  ANALYTICAL 
C  SOLUTION  OF  THE  TATE  EQUATIONS: 

C 

C  dL/dt  -  V  -  U 

C 

C  L  dU/dt  -  -Y/8r 

C  2  2 

C  1/2  er  (U-V)  +  Y  -  1/2  et  V  +  R 

C 

C  dP/dt  -  V 

c 

C  WHERE : 

C  V  -  ROD  VELOCITY 

C  U  -  PENETRATION  VELOCITY 

C  Y  -  TATE  ROD  STRENGTH  PARAMETER 
C  0r  ~  ROD  DENSITY 

C  R  -  TATE  TARGET  STRENGTH  PARAMETER 

C  0t  >  TARGET  DENSITY 

C  L  -  ROD  LENGTH 

C  P  -  PENETRATION 

C  t  -  TIME 

C 


C  ——VARIABLES—— 

C 

C  **»**'DIMENSIONAL  TATE  CONSTANTS**** 

C 

C  R  -  TARGET  RESISTANCE 

C  Y  -  ROD  RESISTANCE 

C  RHOT  -  TARGET  DENSITY 

C  RHOR  -  ROD  DENSITY 

C  Lo  -  INITIAL  ROD  LENGTH 

C  Uo  -  INITIAL  ROD  VELOCITY 

C 

C  *****NON-DIMENSIONAL  CONSTANTS***** 

C 


C  GAMMA 

C  SIGMA 

C  K 

C  A 

C  BIGB 

C  C 

C  M 

C  N 

C  F 

C  ZO 

C  ZX 

c 

C  *****OUTPUT  VARIABLES***** 

C 

C  ZO  -  ARRAY  OF  TRANSFORMATION  VARIABLES,  AT  WHICH  SOLUTION 

C  IS  TO  BE  EVALUATED 

C  TO  -  ARRAY  OF  TIMES,  CORRESPONDING  TO  THE  TRANSFORMATION 
C  Z  VARIABLES 

C  U  -  ROD  VELOCITY 

C  V  -  PENETRATION  VELOCITY 

C  Ldot  -  RATE  OF  ROD  EROSION 

CL-  LENGTH  OF  UNERODED  ROD 

C  PO  -  ARRAY  OF  PENETRATION  VALUES  CORRESPONDING  TO  THE 
C  TRANSFORMATION  Z  VARIABLES 


-  A  CONGLOMERATION  OF  CONSTANTS,  APPEARING  OCCASIONALLY 

-  INITIAL  VALUE  OF  Z  (WHEN  U  -  Uo) 

-  terminal  value  of  Z  (WHEN  V  -  0) 
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XTAIL  -  CCXJBDINATE  OF  THE  REAR  END  OF  THE  ROD 

*****B00KKEEPING  VARIABLES 

MXPS  -  MAXIMUM  NUMBER  OF  TIMES,  LESS  1,  AT  WHICH  SOLUTION  MAY 
BE  EVALUATED 

EPSILON  -  RELATIVE  ERROR  TO  BE  USED  IN  TESTING  FOR  SERIES 
CONVERGENCE  OF  THE  W  INTEGRAL 

NTERMS  -  NUMBER  OF  TERMS  REQUIRED  TO  CONVERGE  THE  W 
INTEGRAL 

NTAVG  -  AVERAGE  NUMBER  OF  TERMS  REQUIRED  TO  CONVERGE  THE  W 
INTEGRAL,  OVER  THE  5  INTEGRAL  EVALUATIONS  REQUIRED. 

NPASS  -  NUMBER  OF  TIMES  AT  WHICH  THE  SOLUTION  IS  EVALUATED 
(MAY  NOT  EXCEED  MXPS+1) 

I  -  A  LOOP  INDEX 

ANSWER  -  A  ONE  CHARACTER  RESPONSE  TO  THE  QUERY  FOR  SCREEN  OR  FILE 
OUTPUT 

*****OTHER  VARIABLES***** 


ZZ  -  SCALAR  VERSION  OF  Z,  CORRESPONDING  TO  THE  SINGLE  Z  ELEMENT 
OF  CURRENT  INTEREST 

DT  -  USER  SPECIFIED  DELTA  TIME,  AT  WHICH  HE  HOPES  TO  HAVE  SOLUTIONS 
EVALUATED  (NOTE  THAT  ACTUAL  SOLUTION  INTERVAL  IS  ONLY 
APPROXIMATELY  THIS  VALUE) 

DZ  -  INCREMENT  IN  Z  CORRESPONDING  ROUGHLY  TO  THE  TIME  INCREMENT  DT 
Ur  -  NON-DIMENSIONAL  ROD  VELOCITY  (U/Uo) 

Vr  -  NON-DIMENSIONAL  PENETRATION  VELOCITY  (V/Vo) 

D  -  AN  EXPONENT  ON  2,  APPEARING  IN  THE  W  INTEGRAL 

TERMxO  -  WHERE  x  IS  1,  2  OR  3...  ARRAYS  CONTAINING  EVALUATIONS  OF 

THE  W  INTEGRAL,  BACH  ARRAY  ELEMENT  WITH  DIFFERENT 
LIMITS  OF  INTEGRATION,  AS  GIVEN  BY  THE  Z  ARRAY 
ROOTGAMMA  -  SQUARE  ROOT  OF  GAMMA,  WHICH  OCCURS  FREQUENTLY 
ROOTTERM  -  A  SQUARE  ROOT  TERM,  APPEARING  OCCASIONALLY 
ROOTZ  -  SQUARE  ROOT  OF  2,  APPEARING  OCCASIONALLY 


C 


integer  MXPS 
PARAMETER  (MXPS>200) 

double  precialon  gairana,  sigma,  K,  A,  BigB,  c,  F,  zO,  zx,  M,  N 

common  /wconstants/  gamma,  sigma,  K,  A,  BigB,  c,  F,  zO,  zx,  M,  N 

double  precision  R,  Y,  rhot,  rhor,  Lo,  Uo 

common  /tconstancs/  R,  Y,  rhot,  rhor,  Lo,  Uo 

double  precision  rootgamma,  z(0:MXPS),  zz,  D,  terml (0:MXPS) , 

&  term2 (0:MXPS) ,  term3 (0 :MXPS) ,  EPSILON,  dt, 

£  dz,  t(0:MXPS),  U,  V,  Ldot,  L,  P(0:MXPS),  xtail, 

&  rootterm,  rootz,  Ur,  Vr,  Cl,  E 

integer  nterms,  ntavg,  npass,  i 
character*!  answer 
PARAMETER  (EPSILON-1 . OD-3) 


1  format 

('  Enter  approximate  dt:  ') 

2  format 

(lx,'  z 

t  U  V 

& 

'dL/dt 

L  P  Xtail') 

3  format 

(f6.2, lp,e9. 

,2,2el0.3,el0.2,2el0.3,ell. 3) 

4  format 

( ' xxxxx' , ' 

8',/, 

& 

'R,  Y,0t,0r 

-  ',lp4el0.3,/. 

& 

'z' 

& 

't',/, 

& 

'U',/, 

& 

'V',/, 

£ 

'dL/dt' ,  /, 

£ 

'L',/, 
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i  'P',/r 

»•.  'Xtall') 

5  format  (Ip/8«10.3) 

6  format  ('  Scraen  Only  or  File  Output  too  (S/F)7  ') 

7  format  Illegal  Choloe. . .Try  again...') 

8  t'ormat  (/, '  Average  nundMr  of  aeriea  terms  used  In  solution  -  ',i3 

a  Number  of  times  at  which  solution  was  evaluated  ■  ',13) 

9  format  <'  Requested  timestep  requires  more  than  ',15, 

a  '  solutions...',/,'  Please  increase  MXPS  parameter') 

10  format  <al) 

11  format  ('  Solutions  have  been  requested  at  more  than',i4,'  times.' 
a  ,/,'  MXPS  parameter  must  be  increased  in  order  to  proceed.') 

12  format  ('  Penetrator  goes  rigid,  since  U  -  V.') 

13  format  ('  Rod  continues  to  erode  in  rigid  target,  since  V  -  0.') 

INPUT  TATE  CONSTANTS  AND  INITIALIZE  CONSTANTS  FOR  CURRENT  MODEL 
call  getconstants 
rootgamna  -  dsqrt (gamma) 

DETERMINE  DESIRED  DT  STEP  SIZE  (NOT  NEEDED  TO  SOLVE  THE  EQUATIONS; 
RATHER,  IT  PROVirES  CONVENIENCE  OF  OUTPUT  SPECIFICATION  TO  THE 
USER) 

write  (*,1) 
read  (*,*)  dt 

INITIALIZE  VARIABLES 
npass  1 

DETERMINE  OUTPUT  MODE  (S  FOR  SCREEN  ONLY;  F  FOR  FILE  OUTPUT  TOO) 

70  write  (*,6) 

read  (*,10)  answer 
CONVERT  ANSWER  TO  UPPER  CASE 
if  (answer  .eq.  's')  then 
answer  ■  'S' 

else  if  (answer  .eq.  '£')  then 
answer  «  'F' 
end  if 

if  (answer  .eq.  'F')  then 

SET  UP  OUTPUT  FILE  'TATE. OUT' 
open  (2, file*' tate.out' , access-' append' ) 
write  (2,4)  R,  Y,  rhot,  rhor 
else  if  (answer  .ne.  'S')  then 
write  (*,7) 
goto  70 
end  if 

FOR  FIRST  ITERATION,  START  AT  ZO 
z(0)  -  zO 

AS  LONG  AS  Z  IS  ABOVE  MINIMUM  LIMIT,  PROCEED  WITH  THE  FOLLOWING: 
100  if  (z(npass'l)  .gt.  zx)  then 

COMPUTE  OZ  REQUIRED  TO  PRODUCE  DESIRED  T  (APPROXIMATE  ONLY) 

BY  EVALUATING  DERIVATIVE  AT  MIDDLE  OF  Z  INCREMENT 
zz  -  z(npass-l) 
do  105  i  -  1,  2 
dz  -  dt 

a  *  F  /  dexp(BigB*zz  -  c/zz) 
a  /  (  zz*M(A-l.)/2.) 

a  +  sigma* (1 . -gamma)  •  zz** ( (A-3 . ) /2 . ) ) 

if  (-dz  .gt.  zz)  then 
zz  -  zz/2. 
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else 

EZ  >  EC  dE/2. 
end  if 

105  continue 
C 

C  DETERMINE  NEW  Z. 

C  IF  Z  FALLS  BELOW  MINIMUM  LIMIT,  RESET  Z  TO  MINIMUM  LIMIT 

z(npa88)  z(npa88-l)  +  ds 
If  (z(npa88)  .It.  zx)  s(npa88)  -  cx 
c 

c  INCREMENT  COUNTER  ON  NUMBER  OF  TIMES  AT  WHICH  SOLUTION  IS 

C  CC»«'UTED 

npass  *  npass  1 

C 

C  CHECK  TO  SEE  THAT  ARRAY  LIMITS  HAVE  NOT  BEEN  EXCEEDED 

if  (npass  .gt.  MXPS)  then 
write  (*,11)  MXPS 
stop 
end  if 
goto  100 
end  if 
C 

C  EVALUATE  T,  FOR  THE  GIVEN  Z,  REQUIRING  SOLUTION  OF  THE 
C  W  INTEGRAL 

D  -  (A-l.)/2. 

call  wintegral  (terml, npass, BigB,c,D,zO,z,ntenns, EPSILON) 
ntavg  -  nterms 
C 

D  -  (A-3.)/2. 

call  wintegral  (tenn2, npass, BigB, c,D, zO, z, nterms, EPSILON) 
ntavg  «  ntavg  +  nterms 
C 

do  110  i  -  0,  (npass-l) 

t(i)  -  (terml(i)  +  sigma* (1 . -gamma) *term2 (i) )  /  F 
110  continue 
C 

C  EVALUATE  PENETRATION,  REQUIRING  SOLUTION  OF  W  INTEGRAL 
D  -  (A  ) /2. 

call  wintegral  (terml, npass, BigB, c,0, zO, z, nterms, EPSILON) 
ntavg  ■  ntavg  +  nterms 
C 

D  -  (A-2. ) /2. 

call  wintegral  (term2, npass, BigB, c,0, zO, z, nterms, EPSILON) 
ntavg  -  ntavg  +  nterms 

C 

D  -  (A-4 . ) /2. 

call  wintegral  (term3, npass, BigB, c, D, zO, z, nterms, EPSILON) 
ntavg  ■  (ntavg  +  nterms)  /  5 
C 

do  115  i  ■  0,  (npass-l) 

P(i)  -  Uo  /  F  *  ( 

&  (1-rootgamma) / (2 . *rootgan»na* (1 . -gamma) )  *  terml (i) 

&  sigma  *  term2(i) 

&  sigma*  *2* (1 . +rootgamma) * (1 . -gamma) / (2 . *rootgamma) 

&  *  term3 (i) ) 

115  continue 
C 

C  SET  UP  HEADER  COLUMNS 

write  (*,2) 

C 

C  DETERMINE  REST  OF  PARAMETERS 

do  120  i  “  0,  (npass-l) 
zz  -  z(i) 
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EVALUATE  U,  FROM  THE  GIVEN  2 

Ur  *  (zz  -  sigma*  (1  .-gansna) )  /  <2.*tootgainma*dsqrt  (zz) ) 

U  -  Ur  *  Uo 

EVALUATE  V  PROM  THE  C<»1PUTB0  U 

Vr  •  (Ur  -  dsqrt  (gainma*Ur**2  +  aigma*  <1 . -gamma) ) )  /  (1. -gamma) 
V  -  Vr  *  Uo 

EVALUATE  Ldot  FROM  COMPUTED  U  AMD  V 
Ldot  -  V  -  U 

EVALUATE  ROD  LENGTH  FROM  COMPUTED  U 

rootterm  -  dsqrt (gamma*Ur**2  +  sigma* (1. -gamma) ) 

rootz  Ur*rootgaitina  +  rootterm 

L  ••  Lo  •  dexp (Ur*  (rootterm  -  Ur*gamma)  /  (2 . *K*  (1  .-gamma) ) ) 
i  *  (rootz) **A  /  (dexp(N)  *  M) 


COMPUTE  TAIL  LOCATION  OF  ROD 
xtail  -  P(i)  -  L 

WRITE  TO  SCREEN 

write  (*,3)  z(i),  t(i),  U,  V,  Ldot,  L,  P(i),  xtail 
OUTPUT  RESULTS 
if  (answer  .eq.  'F')  then 

WRITE  TO  FILE 

write  (2,5)  z(i),  t(i),  U,  V,  Ldot,  L,  P(l),  xtail 
end  if 
120  continue 

COMPLETE  THE  PROBLEM  TO  ALLOW... 
i  -  npass 
z(i)  -  0. 

U  -  0. 

V  -  0. 

Ldot  -  0. 

. . .RIGID  BODY  ROD  TO  COME  TO  STOP,  IF  R  <  Y 
if  (sigma  .It.  0.)  then 
write  (*,12) 

L  -  L 

P(i)  -  P(i-l)  +  L/gamma  *  dlog(y/R> 

Cl  -  dsqrt  (  2.*K'*R/(gamma*Yi  ) 
t(i)  "  t(i-l)  +  2 . *L/ (gamma-Uo*Cl)  * 

&  datan  (Ur/Cl) 

...OR  NON-PENETRATING  ERODING  ROD  TO  COMB  TO  STOP,  IF  R  >  Y 

else 

write  (*,13) 

L  -  L  *  dexp  ((Y-R)/Y) 

P(i)  -  P(i-l) 
z (0)  -  0 . 

BigB  -  l./(2,*K) 
c  -  0, 

E  -  2. 

call  wintegrate  (terml, 1 , BigB, c, E, Ur, z, nterms, EPSILON) 
t(i)  »  t(i-l)  -  L/(2.*K*Uo)  *  terml (0) 
end  if 

xtail  -  P(i)  -  L 
C  WRITE  TO  SCREEN 

write  (*,3)  z(i),  t(i),  U,  V,  Ldot,  L,  P(i),  xtail 
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C  OUTPUT  RBSULTS 

if  (answer  .eq.  'P')  then 

WRITE  TO  FILE 

write  (2,5)  *(i),  t(i),  U,  V,  Ldot,  L,  P(i),  *tail 
end  if 

FINISHED  ITERATING...  CLOSE  UP  SHOP, 
write  (*,8)  ntavg,  npass 
if  (answer  .eq.  'F')  close  (2) 
stop 
end 

********************************** ************************************* 

subroutine  getconstants 

ACQUIRE  THE  TATE  CONSTANTS  NEEDED  BY  THE  CODE,  AND  INITIALIZE 
CORRESPONDING  CONSTANTS  FOR  PRESENT  SOLUTION 

—.—VARIABLES—— 

*****DIMENSI0NAL  TATE  CONSTANTS***** 

R  -  TARGET  RESISTANCE 
y  -  ROD  RESISTANCE 
RHOT  -  TARGET  DENSITY 
RHOR  -  ROD  DENSITY 
LO  -  INITIAL  ROD  LENGTH 
Uo  -  INITIAL  BOD  VELOCITY 

•****NON-DIMENSIONAL  CONSTANTS***** 


GAMMA 

SIGMA 

K 

A 

BIGB 

C 

M 

N 

ZO  -  INITIAL  VALUE  OF  Z  (WHEN  U  -  Uo) 
ZX  -  TERMINAL  VALUE  OF  Z  (WHEN  V  -  0) 

♦*»* ‘OTHER  VARIABLES***** 


•.  -  MINIMUM  VALUE  OF  U,  WHEN  V  -  0 


Di'ble  precision  gamma,  sigma,  K,  A,  BigB,  c,  F,  zO,  zx,  M,  N 
Cs^mmon  /wconstants/  gamma,  sigma,  K,  A,  BigB,  s,  F,  zO,  zx,  M,  N 
double  precision  R,  Y,  rhot,  rhor,  Lo,  Uo 
common  /tconstants/  R,  Y,  rhot,  rhor,  Lo,  Uo 
double  precision  Umin 


1  format 

2  format 

3  format 

4  format 
6 

5  format 


('  Enter  the  target  4  rod  resistances  (R  «  Y) :  ') 
('  Enter  the  target  i  rod  densities:  ') 

('  Enter  the  initial  rod  length  and  velocity:  ') 
dp,'  The  target  resistance,' ,el0. 3, 

',  must  not  equal  the  rod  resistance.') 
dp,'  T.he  initial  rod  velocity, '  ,el0 . 3, 

',  must  exceed  the  minimum  value, ' , elO . 3, ' . ' ) 


READ  PARAMETERS  REQUIRED  IN  TATE  EQUATIONS 
write  (*,1) 
read  (*,*)  R,  Y 
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writ*  i*,2) 

r**d  (*,*)  rhot»  rhor 

writ*  (*,3) 

r*ad  (*,*)  Lo,  Uo 

CHECK  FOR  DATA  INCONSISTBMCIBS 

TARGET  RESISTANCE  MUST  NOT  EQUAL  ROD  RESISTANCE,  FOR 
SOLUTION  TO  TATE  EQUATIONS  TO  BE  MEANINGFUL 
if  (R  .eq.  Y)  th4*n 
writ*  (*,4)  R 
atop 
•nd  if 

Umin  OCCURS  WHEN  PENETRATION  VELOCITY  V  •  0  FOR  R  >  Y,  AND  U  -  V 
FOR  R  <  Y.  IT  ROUGHLY  CORRESPONDS  TO  THAT  VELOCITY  WHERE  THE  ROD 
BOUNCES  OFF  THE  TARGET,  WITHOUT  MAKING  PENETRATION.  MAKE  SURE 
THAT  UO  EXCEEDS  THE  VALUE  OF  Umin. 
sigma  -  2.  *  (R-Y)  /  (rhor*Uo**2) 
gamma  -  rhot  /  rhor 
if  (sigma  .gt.  0.)  than 
Umin  -  Uo  *  dsqrt (sigma) 

*ls* 

Umin  >  Uo  *  dsqrt  (-signta/gamna) 

«nd  if 

if  (Uo  .It.  Umin)  than 
write  (*,S)  Uo,  Umin 
stop 
end  if 

INITIALIZE  REST  OF  CONSTANTS  (DERIVED  FROM  TATE  CONSTANTS) 

call  initialize 

return 

end 

subroutine  initialize 

INITIALIZE  CONSTANTS  USED  IN  SOLUTION  OF  TATE  PROBLEM, 

GIVEN  THE  TATE  CONSTANTS 

.——VARIABLES—— 

*****DIMENSI0NAL  TATE  CONSTANTS***** 

R  -  TARGET  RESISTANCE 
Y  -  ROD  RESISTANCE 
RHOT  -  TARGET  DENSITY 
RHOR  >  ROD  DENSITY 
Lo  >  INITIAL  ROD  LENGTH 
Uo  -  INITIAL  ROD  VELOCITY 

*  *  *  *  *NON-DIMENS lONAL  CONSTANTS*  *  *  *  * 


GAMMA 

SIGMA 

K 

A 

BIGB 

C 

M 

N 
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ZO  -  INITIAL  VALUE  OF  Z  (WHEN  U  •  Uo) 

ZX  -  TERMINAL  VALUE  OF  Z  (NHEN  V  -  0) 

i»****OTHER  VARIABLES***** 

ROOTGAMMA  -  SQUARE  ROOT  OF  GAMMA,  WHICH  OCCURS  FREQUENTLY 
ROOTTERM  -  A  SQUARE  ROOT  TERM,  APPEARING  OCCASIONALLY 
ROOTZO  -  SQUARE  ROOT  OF  ZO,  APPEARING  OCCASIONALLY 
UdOCOr  -  ORIGINAL  ACCELERATION  OF  ROD  /  Uo 

double  precision  gaaam,  eigme,  K,  A,  BlgB,  c,  T,  zO,  ^x,  M,  N 
conmon  /wconetante/  gemne,  eigne,  K,  A,  BigB,  c,  F,  zO,  z'^,  M,  N 
double  precision  R,  Y,  rhot,  rhor,  Lo,  Uo 
conroon  /tconstants/  R,  Y,  rhot,  rhor,  Lo,  Uo 
double  precision  rootganma,  roottem,  rootzO,  Udotor 

GAMMA  DEFINED  IN  ROUTINE  GBTC<»«STANT8. . .NO  NEED  TO  REINITIALIZE 
ganroa  ~  rhot  /  rhor 

TEMPORARY  VARIABLE... 
rootgamca  -  dsq't (gamna) 

SIGMA  defined  in  ROUTINE  OETCONSTANTS. . .NO  NEED  TO  REINITIALIZE 
Sigma  -  2.  *  (R-Y)  /  (rhor*Uo**2) 

K  -  Y  /  (rhor*Oo**2) 

A  -  sigma  /  <2.  *  K  *  rootgansna) 

BigB  -  1.  /  (8.  *  K  *  rootgamna  *  (rootganma  t  1.)) 

c  -  sigma**2  *  (1,  -  gamma)  •  (rootganma  +1.)  / 

i  (6.  *  K  *  rootganma) 

ZO  IS  THE  INITIAL  VALUE  OF  Z,  CORRESPONDING  TO  THE  SITUATION 
OF  U  -  Uo. 

rootterm  -  dsqrt (gamma  t  sigma* (1 . -gamma) ) 
rootzO  *  rootganma  t  rootterm 
rO  -  rootz0**2 

M  AND  N  ARE  CONSTANTS  APPEARING  IN  THE  EQUATION  FOR  dU/dt  (AND  THUS 
THE  EQUATION  FOR  U(t) . 

M  ••  (rootzO)  **A 

N  ■  (rootterm  -  gamma)  /  (2.*K* (l.-gamma) ) 

Udotor  “  -Y  /  (rhor  *  Lo  *  Uo) 

F  -  4.  *  Udotor  *  M  *  rootgamma  *  dexp(N  -  sigina/(4.  *  K) ) 
if  (R  .gt.  Y)  then 

( 

FOR  R  >  Y: 

ZX  IS  THE  TERMINAL  VALUE  OF  Z,  WHEN  V  -  0.  FROM  'niE  EQUATION, 

WE  SEE  THAT  V  -  0  WHEN  U/Uo  -  ROOT (SIGMA).  WHEN  U  IS  THIS  VALUE, 
Z  TAKES  ON  THE  FOLLOWING  VALUE: 

ZX  ■  sigma  *  (rootgannna  +  l.)**2 
else 

FOR  R  <  Y; 

ZX  IS  THE  TERMINAL  VALUE  OF  Z,  WHEN  V  -  U.  FROM  TATE  EQUATION, 

WE  SEE  THAT  V  -  U  WHEN  U/Uo  -  ROOT (-SIGMA/GAMMA) .  WHEN  U  IS 
THIS  VALUE,  Z  TAKES  ON  THE  FOLLOWING  VALUE: 

ZX  -  -sigma  *  (rootgamma  +  l.)**2 
end  if 
C 

return 

end 

c*********************** ********************************* ****•***•**•*•• 

subroutine  wintegral  (w,  npass,  B,  C,  D,  zl,  z2,  n,  eps) 
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100 


SUBROUTINE  TO  IVM.UATB  THB  N  XNTBQRAL: 
b2 

(  D  -1 

t  z  exp  (Bz  -  Cz  )  dz 

) 

zl 

TO  DO  SO,  TRANSFORM  THB  INTBORAL  NXTH  THB  FOLLOMXNO  TRANSFORMATION 
AND  INTEGRATE  THAT  FUNCTION,  AS  FOLLOWS: 

D+1  1/(D+1)  D 

let  y  *  z  THBN  z  *  y  AND  dy  •  (O'*-!)  z  dz. 

LET  E  -  1/(D+1).  THB  INTEGRAL  THEN  BECOMES 

y2 

(  E  -E 

E  I  exp  (By  -  cy  )  dy 
) 

yi 


...i»VARIABLBS— 

**** •bookkeeping  variables 

MXPS  -  MAXIMUM  NUMBER  OF  TIMES,  LESS  1,  AT  WHICH  SOLUTION  CAN 
BE  EVALUATED 

EPS  -  RELATIVE  ERROR  TO  BE  USED  IN  TESTING  FOR  SERIES 
CONVERGENCE  OF  THB  N  INTEGRAL 
N  >  NUMBER  OP  TERMS  REQUIRED  TO  CONVERGE  THB  W 
INTEGRAL 

I  -  A  LOOP  INDEX 

NPA3S  -  NUMBER  OF  Z  VALUES  AT  WHICH  TO  EVALUATE  SOLUTION 

** •••other  variables***** 


MO 


Zl 

Z2() 

B 

c 

0 
E 
Yl 
Y2  0 
DPLUSl 


ARRAYS  CONTAINING  EVALUATIONS  OF  THB  M  INTEGRAL,  BACH 
ARRAY  ELEMENT  WITH  DIFFERENT  LIMITS  OF  INTEGRATION,  AS 
GIVEN  BY  THE  Z2  ARRAY 
LOWER  LIMIT  OF  INTEGRATION 

ARRAY  CONTAINING  UPPER  LIMITS  OF  INTEGRATION 
CONSTANT  MULTIPLIER  OF  Z  IN  EXPONENTIAL  TERM  OF  W 
INTEGRAL 

CONSTANT  MULTIPLIER  OF  1/Z  IN  EXPONENTIAL  TERM  OF  W 
INTEGRAL 

EXPONENT  OF  Z  IN  N  INTEGRAL 
TRANSFORMED  VARIABLE  EQUAL  TO  l/(D-t-l) 

TRANSFORMED  Zl  LIMIT  OF  INTEGRATION 
TRANSFORMED  Z2()  ARRAY  LIMITS  OF  INTEGRATION 
INTERMEDIATE  VARIABLE 


integer  MXPS 
PARAMETER  (MXPS-200) 
integer  n,  i,  npasa 

double  precision  w(0:MXPS),  B,  C,  D,  zl,  z2(0:HXPS),  eps 
double  precision  B,  yl,  y2(0:HXPS),  Dplusl 
Oplusl  «  D+1. 
yl  ■  zl**Dplusl 
do  100  i  -  0,  (npass'l) 
y2(i)  -  z2(i) **Dplusl 
E  -  1. /Dplusl 
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call  wlntagrat*  (w,  npaea,  B,  C,  t,  yl,  y2,  n,  ape) 
raturn 

and  ^  ^ 

subroutine  wintagrata  (w»  npaaa,  B,  C,  E,  yl/  y2,  n/  aps) 

FUNCTION  TO  EVALUATE  THE  H  INTEGRAL  <HITH 
TRANSFORMATION)/  BY  SERIES  EXPANSION: 


B  I  exp  (By  -  Cy  )  dy 
) 

yl 

THE  SERIES  EXPANSION  OF  THE  EXPONENTIAL  13: 

2  3 

exp  (x)  -  1  +  X  +  X  /2!  +  X  /3!  ... 


Ill  t 

TERM  *  0  1  2  3 

WHEN  THE  ABSOLUTE  VALUE  OF  A  GIVEN  TERM  OF  THIS  EXPANSION  IS 
LESS  THAN  THE  ERROR  "apa**  TIMES  THE  CURRENT  FUNCTION  VALUE/ 

THE  SERIES  IS  TRUNCATED.  THE  NUMBER  OF  TERMS  REQUIRED  TO 
SUM  WITHIN  THIS  ERROR  VALUE  IS  RETURNED  IN  VARIABLE  "n". 

NOTE  THAT  INTEGRAL  OF  TERM  0  OF  THE  EXPANSION  IS  (y2-yl)  WHEN 
EVALUATED  BETWEEN  yl  AND  y2. 

REMEMBER  ALSO/  THAT  TERM  "x«  IN  THIS  EXPANSION  IS  REALLY  A  TERM 
LIKE  (A  <>■  B)  .  THUS/  THE  POLYNOMIAL  EXPANSION  OF  TERMS  IN  THE 
EXPONENTIAL  ARB: 

n  n  n  k  (n<*k) 

(A  4  B)  -  SUM  (  )  A  B 

k-0  k 

WHERE: 

n  n! 

(  )  - . . 

k  (n-k) !  k! 

integer  MXPS 
PARAMETER  (MXPS-200) 
integer  n,  k/  npasa/  i 

double  precision  m(0:MXPS)/  B/  C,  Z,  yl#  y2(0:MXPS)/  epa 
double  precision  f,  fi,  exponent,  expz, 
t  value (0 :tnPS) ,  tvalue (0:MXPS) ,  minlimit 

W  INTEGRAL  ONLY  INTENDED  FOR  POSITIVE  LIMIT  OF  INTEGRATION 
minliinit  ■  yl 
do  80  i  -  0/  (npass-l) 

80  if  (y2(i)  .It.  minlimit)  minlimit  •  y2(i) 

C  CHECK  TO  SEE  IF  INTEGRATION  LIMITS  ARE  OK 
if  (minlimit  .It.  0.)  then 

write  (*/'("  Negative  integration  limits  not  allowed.")') 
stop 
end  if 
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c 


INITIXLIZB  COUNTER  FOR  #  OF  TERMS  OF  EXPONENTIAL  EXPANSION  USED 
n  «  0 


EVALUATE  THE  INTEGRAL  OF  THE  FIRST  TERM  DIRECTLY: 
do  90  i  -  0,  (npass-l) 

90  N(i)  -  (y2<i)  -  yl) 

FOR  EVERY  ADDITIONAL  TERM  OF  THE  EXPONENTIAL  EXPANSION,  DO  THE 
FOLLOWING: 

INCREMENT  TERM  COUNTER 
100  n  ■  n  +  1 

INITIALIZE  INTEGRAL  ASSOCIATED  WITH  THIS  EXPONENTIAL  TERM  TO  ZERO 
do  110  1-0,  (npa88-l) 

110  tvalue(l)  -  0. 


LOOP  OVER  BACH  POLYNOMIAL  EXPANSION  TERM,  FOR  THIS  EXPONENTIAL 
TERM 

do  200  Ic  -  0,  n 

GET  THE  CONSTANT  ASSOCIATED  M/  THIS  EXPONENTIAL/POLYNOMIAL 
PRODUCT  TERM... IT  IS  GIVEN  BY: 
f  -  B^Jc  *  (-C)''(n-k)  /  (n-k)  !  k! 

THE  FOLLOWING  TECHNIQUE  REDUCES  PROBABILITY  OF  NUMERICAL 
OVERFLOW,  BY  INTERSPERSING  DIVISION  AND  MULTIPLICATION 
f  -  1. 

do  120  i  -  1,  maxO (k, (n-k) ) 

£i  -  dfloatd) 

if  (i  .18.  k  )  f  -  £  *  b  /£i 
if  (i  .18.  n-k)  £  -  £  *  (-c)/£i 
120  continue 

if  (f  .n8.  0.)  th8n 

GET  THE  EXPONENT  ON  Y  ASSOCIATED  W/  THIS  POLYNOMIAL  EXPANSION 
TERM 

8xpon8nt  -  E  *  (2*k  -  n) 


BEGIN  THE  INTEGRATION  OF  THIS  POLYNOMIAL  EXPANSION  TERM: 


C 

C 


y2 

{  n  k  (n-k)  E(2k-n)  y2 

I  (-1)  B  C  y  ( 

I  - dy  .  I 

)  (n-k) !  k!  ) 

yl  yl 

if  (exponent  .ne.  -1.)  then 

IF  EXPONENT  NOT  EQUAL  TO  -1,  THEN 


exponent 
f  y  dy 


y2 

(  exponent  £  (exponent+1)  Iy2 

I  fy  dy- -  y  I 

)  (exponent  +1)  lyl 

yl 


FOR  NORMAL  POLYNOMIAL  INTEGRATION... 

GET  THE  EXPONENT  AFTER  INTEGRATION 
exp2  -  exponent  +  1 . 

GET  THE  TOTAL  COMBINED  CONSTANT  AFTER  INTEGRATION 

f  -  £  /  exp2 

do  150  i  -  0,  (npa33-l) 

EVALUATE  THE  INTEGRAL  BETWEEN  Yl  AND  Y2 
value (i)  -  f  •  (y2(i)**exp2  -  yl**exp2) 
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LUMP  TMZ8  IMTSOltATION  VALUt  FOR  THIS  POLYNOHZAL  TBRM  INTO 
THE  SUN  FOR  TRX  KXP(»ntNTZAL  TERN 
tvalu«(i)  •  tvalued)  *  valued) 

150  continue 

elae 

OTHERNISB,  IF  EXPONENT  -  -1,  THEN 

y2 

(  e)q>onent 

I  f  y  dy  •  f  ln{y2/yl) 

) 

yi 

FOR  LOGARITHMIC  INTEGRATION... 

THE  TOTAL  COMBINED  CONSTANT  AFTER  INTEGRATION  IS  JUST  F 
do  160  i  -  0,  (npaes-l) 

EVALUATE  THE  INTEGRAL  BETWEEN  Yl  AND  ¥2 
valued)  ■  f  *  dlog{y2 d) /yl) 

LUMP  THIS  INTEGRATION  VALUE  FOR  THIS  POLYNOMIAL  TERM  INTO 
THE  SUM  FOR  THE  EXPONENTIAL  TERM 
tvalued)  ~  tvalued)  *  valued) 

160  continue 

end  if 
end  if 

200  continue 

LUMP  THIS  INTEGRATION  VALUE  FOR  THE  EXPONENTIAL  TERM  INTO 
THE  SUM  FOR  THE  N  INTEGRAL 
do  250  i  “  0,  (npass-1) 

250  w(i)  -  w(i)  tvalue(i) 

CHECK  FOR  SERIES  CONVERGENCE  AGAINST  USER  SPECIFIED  EPSILON; 

IF  TERM  DOESN'T  CONVERGE  FOR  ANY  OF  THE  LIMITS  OF  INTEGRATION, 
THEN  COMPUTE  ANOTHER  EXPONENTIAL  SERIES  TERM  FOR  *ALL*  THE 
LIMITS  OF  INTEGRATION 
do  300  i  •  0,  (npas8-l) 

300  if  (dabs  (tvalue  d) )  .gt.  dabs  (eps*wd) ) )  goto  100 

MODIFY  INTEGRAL  VALUE  TO  ACCOUNT  FOR  E  CONSTANT  IN  FRONT 
OF  INTEGRAL 

do  350  i  >  0,  (npass-l) 

350  w(i)  -  E  *  w(i) 
return 
end 
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1  Cenual  Intelligence  Agency 
Office  of  Central  Reference 
Dissemination  Branch 
Room  GE47  HQS 
Washington,  DC  20502 

1  Advanced  Technology,  Inc. 

ATTN;  John  Adams 
PO  Box  125 

Dahlgren,  VA  22448-0125 


No.  of 

Copies  Oraaniatton 

14  Sandia  National  Laboratories 
ATTN:  Division  9122, 

Robert  O.  Nellums 
Jim  HicloHson 
Vince  luk 
Division  1533, 

Marlin  Kipp 
Alien  Robinson 
Division  2512, 

WilliaiTl  J.  Andrzejewski 
Don  Marchi 
Division  1551, 

R.  Graham 
R.  Lafarge 
M.  VigU 
R.  Sandoval 
J.  Asay 
R.  Longcope 
M.  Porrestal 
POBox  5800 
Albuquerque,  NM  87185 

11  Director 

Los  Alamos  National  Laboratory 
ATTN:  MSK574, 

O.  E.  Cort 
Tony  RoUeu 
Mike  Burkett 
MS  P940, 

Robert  Karpp 
MS  K557  N-6, 

Rudy  Henninger 
MS  G740. 

Roy  Greiner 
B214  T.14, 

James  P.  Ritchie 
MS  G787. 

John  Bolsiad 
J.  Walsh 
C.  Mautz 
Tech  Library 
PO  Box  1663 
Los  Alamos,  NM  87545 

I  Explosive  Technology 

ATTN:  Michael  L.  Knacbel 
POBox  KK 
Fairfield,  CA  94533 
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No.  of 

No.  of 

CPDiss 

OoBwi^tias 

Cooies 

OntanizatiQB 

14 

Lawrence  Livennore  National  Laboratories 
ATIN:  L-122. 

Bany  R.  Bowman 

Ward  Dixon 

Raymond  Pierce 

1 

Rockwell  Missile  Systems  Division 
ATTN:  Terry  Neuhart 

1800  Satellite  Blvd. 

Duluth.  GA  30136 

Rusaell  Rosinsky 

Owen  J.  Alford 

1 

Rockwell  IntlVRocketdyne  Division 
ATTN:  James  Moldenhauer 

Diana  Stewart 

Tony  Vidlak 

L-290. 

6633  Canoga  Ave  (HB  23) 

Canoga  Park,  CA  91303 

Albert  Holt 

John  E.  Reaugh 

L-352. 

David  Wood 

L-874. 

Robert  M.  fCuklo 

2 

McDonnell-Douglas  Helicopter 

ATTN:  Loren  R.  Bird 

Lawrence  A.  Mason 

5000  E.  McDoweU  Rd.  (MS  543-D216) 
Mesa,AZ  85205 

MS-35. 

Thomas  McAbee 

Michael  J.  Murphy 

Technical  Library 

PO  Box  808 

1 

S-Cubed 

ATTN:  Dr.  R.  T.  Sedgwick 

PO  Box  1620 

U  Jolla.  CA  92083-1620 

Livermore,  CA  94550 

2 

Orlando  Technology,  Inc. 

ATTN:  Dan  Matuska 

2 

Southwest  Research  Institute 

ATTN:  C.  Anderson 

A.  Wenzel 

PO  Drawer  28255 

J.  Osbom 

PO  Box  855 

Shalimar,  FL  32579 

San  Antonio,  TX  78228-0255 

3 

Kaman  Sciences  Corporation 

ATTN:  D.  Barnette 

3 

Battelle  •  Columbus  Laboratories 

ATTN:  R.  Jameson 

S.  Golaski 

Technical  Library 

505  King  Avc. 

D.  Elder 

P.  Russell 

PO  Box  7463 

Colorado  Springs,  CO  80933-7463 

Columbus.  OH  43201 

2 

Physics  International 

ATTN:  Ron  Funston 

1 

Defense  Technology  International.  Inc. 

ATTN:  D.  E.  Ayer 

The  Stark  House 

22  Concord  Street 

Nashua,  NH  03060 

2 

Lamont  Garnett 

2700  Merced  St. 

PO  Box  5010 

San  Leandro,  CA  94577 

Lockheed  Missile  &  Space  Co.,  Inc. 

4 

California  Research  and  Tech.  Corp. 

ATTN:  Roland  Franzen 

Dennis  Orphal 

Ron  E.  Brown 

Mark  Majerus 

5117  Johnson  Drive 

Pleasanton.  CA  94566 

ATTN:  S.  Kusumi  (0-81-11.  Bldg  157) 
Jack  Philips  (0-54-50) 

PO  Box  3504 

Sunnyvale.  CA  94088 
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No.  of 

Copies  Organizaiion 

1  Lockheed  Missile  &.  Space  Co..  Inc. 
ATIN:  Richard  A.  Hoffman 
Santa  Cniz  Facility 

Empire  Oiade  Ro^ 

Santa  Cniz,  CA  9S060 

2  Mason  &  Hanger  •  Silas  Mason  Co. 
ATTN:  Thomas  J.  Rowan 

Christopher  Vogt 
Iowa  Army  Ammunition  Plant 
Middletown.  lA  S2638-9701 

1  Lockheed  Engineering  &  Space  Sciences 
ATTN:  Ed  Cykowaki.  MS  B-22 

2400  NASA  Road  1 
Houston,  TX 

2  Dyiia  East  Corporation 
ATTN:  P.  C.  Chou 

R.  Ciccarelli 
3201  Arch  St. 

Philadelphia,  PA  19104 

2  E.  I.  DuPont  De  Nemoura  &  Company 
ATTN:  B.  Scott 

L.  Mina 

Security  Director,  Legal  Department 
PO  Bos  1635 
Wilmington,  DE  19899 

1  Aerojet  Electro  Systems  Company 

ATTN:  Warhead  Systems, 

Df.  J.  Caileone 
PO  Box  296 
. ‘zusa,  CA  91702 

1  Physics  International  Company 

Tactical  Systems  Group 
Eastern  Division 
PO  Box  1004 

Wadsworth,  OH  44281-0904 

3  Alliant  Techsystems,  Inc. 

Defense  Systems  Division 
ATTN:  G.  Johnson 

J.  Houlton 
N.  Berkliolu 
7225  Northland  Drive 
Brooklyn  Park.  MN  55428 


No.  of 

Copies  Organization 

1  Martin  Marietta  Aerospace 

ATIN:  Donald  R.  Bragg 
PO  Box  5837,  MP  109 
Orlando,  FL  32855 

1  SRI  International 

ATTN:  Dr.  L.  Seaman 
333  Ravenswnod  Ave. 

Menlo  Park,  CA  94025 

1  Northrop  Corpoation 

Electro-Mechanical  Division 
ATIN:  Donald  L.  Hall 
500  East  Orangethorpe  Ave. 

Anaheim,  CA  92801 

3  Boeing  Aerospace  Co. 

Shock  Physics  &  Applied  Math 
Engineering  Technology 
ATTN:  R.  Helzer 
J.  Shrader 
T.  Murray 
PO  Box  3999 
Seattle.  WA  98124 

1  McDonnell  Douglas  Astronautics  Company 

ATTN:  Bruce  L.  Coopa 
5301  Bolaa  Ave 
Huntington  Beach,  CA  92647 

1  D.R.  Kennedy  and  Associates  Inc. 

ATTN:  Doriald  Kennedy 
PO  Box  4003 

Mountain  View,  CA  94040 

1  University  of  Colorado 

ATIN:  Timothy  Maclay 
Campus  Box  431.  NNT  341 
Boulder,  CO  80309 

1  New  Mexico  Institute  Mining  &  Tech. 

ATTN:  David  J.  Chavez 
Campus  Station,  TERA  Group 
Socorro,  NM  87801 
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No.  or 


Copies 

1  UniverslQr  of  Dayton 
Reseatch  Inaitote 

Document  Control,  Wright  Brothers  Bnnch 
ATIN:  Dr.  S.  J.  Bless 
PO  Box  283 
Dayton.  OH  43409 

4  Univerhty  of  Delaware 

DqK.  of  Mechanical  Engineering 
ATTN:  Prof.  J.  Vinson 
Prof.  D.  Wilkins 
Prof.  J.  Oilkqiie 
Dean  R.  B.  P^ 

Newark,  DE  19716 
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£SBiSI  QrtMrfmdni}  Cooiet 

2  Defense  Reses'ch  EitsbliAffleot  SofBsld  2 

ATIN:  Chris  Weickeit 
DevkIMsckty 

Ralston.  Albeita.  TQI 2NO  Ralston 
CANADA 

1  Defense  Resevch  EstaUlshineot  Valcanier 

ATIN:  Noiheri  Gass  1 

PO  Box  8800 

Couicelette,  PQ,  GOA  IRO 
CANADA 

1  Canadiui  Arsenals.  LTD 
ATTN;  Piene  Pelletier 

5  Moniee  des  Arseoaux  1 

VUUe  de  Gardeur,  PQ.  J5Z2 

CANADA 

2  Ernst  Mach  Institute 
ATIN:  A.J.  StUp 

V.  Kohler 

Eckmtrasae4  4 

D-7800  Freiburg  i.  Br. 

GERMANY 

3  lABG 

ATIN:  H.  J.  Raatschen 

W.  Schittke 
F.  Scharppf 

Einsteinstrasse  20 

D-8012  Ouobrun  B.  Muenchen 

GERMANY  1 

I  Royal  Aimameni  RAD  Establishment 
ATTN:  Ian  CulUi 
Fort  Halstead 

Sevenoaks,  Kent  TN14  7BJ 
ENGLAND 

I  Centre  d’Etudes  de  Oramat 
ATTN:  Gerald  Solve 
46500  Gramat 
FRANCE 

1  PRB  S.A. 

ATTN:  M.  Vansnick 
Avenue  de  Tervueren  168,  Bie.  7 
Brussels,  B*11S0 
BELGIUM 


OaanlKation 

AB  Bofon/Anunonition  Division 
ATTN  :  JanHassUd 

Anders  Noidell 
BOX  900 
S<691  80  Bofon 
SWEDEN 

Messerschmlu-BOOcow'Blohin 
Dynamics  Division 
ATTN;  Manfred  Held 
PO  Box  1340 
D^98  Schrobenhausen 
GERMANY 

TNOPrins  Maurits  Laboratory 
ATTN:  B.  J.  Pasman 
POBOX4S 
2280  AA  Rii)swijk 
Lange  Kleh^  137 
*»B  NETHERLANDS 

IS 

Al  u.  Weihrauch 
H.  J.  Ernst 
P.  Y.  Chanterei 
V.  Schlr., 

F68301  Saint'Louis 
Cddex,  12.  me  de  ITndusirie 
BP.  301 
FRANCE 

IFAM 

ATTN;  Lodiar  Meyer 
Lesumer  Heerstrasse  36 
Bremen  77 
GERMANY 
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USER  EVALUATION  SHEBT/CHANGE  OF  ADDRESS 


This  Utboritofy  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes. 
Your  commenta/uiswets  to  die  ItemiAlu^oas  below  will  aid  us  in  our  efforts. 

1.  BRL  Report  Number  BRL-TR-3180 _ ntre  nf  Bupnit  pynpiKRR  1990 

2.  Date  Report  Received _ 

3.  Does  ddt  report  satisfy  a  need?  (Oonunent  on  purpose,  related  project,  or  other  area  of  interest 

for  which  the  rqioit  win  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source 
of  ideas,  etc.) _ 


S,  Has  the  information  in  this  report  led  to  any  qusniitaiive  savings  as  far  as  man-hours  or  dollars 
saved,  operating  costs  avoided,  or  efficiencies  achieved,  etc?  if  so,  please  elaborate. _ 


6.  General  Commenu.  What  do  you  think  should  be  changed  to  improve  future  nepoas?  (Indicate 
changes  to  organization,  technical  content,  fonnat.  etc.)  _ _ 


Name 


CURRENT  Organization 

ADDRESS  _ 

Address 


City,  Slate,  Zip  Code 

7,  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  New  or  Corrcti 
Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below 


Name 


OLD  Organization 

ADDRESS  _ _ 

Addres-s 


City,  Stale,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  staple  or  tape  dosed,  and  mail.) 
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OFFICIAL  BUSINESS 


NOPOiTAOE 
NECEMMY 
IF  MAUD 
MTHE 

UMTEOrATEt 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  No  0001,  APG.  MO 


POSTAGE  WUl  K  PAID  BY  ADDAESSEE 


Director 

L'.S.  Army  Ballistic  RcM.-rch  l  iilmratory 

ATTN:  SLCBR-DD  T 

Atxirdccn  Proving  Ground.  MD  21005-‘J989 


FOLD  lILKi: 


